home *** CD-ROM | disk | FTP | other *** search
/ Die Ultimative Software-P…i Collection 1996 & 1997 / Die Ultimative Software-Pakete CD-ROM fur Atari Collection 1996 & 1997.iso / g / gnu_c / pmlsrc23.zoo / pmlsrc / acos.c next >
Encoding:
C/C++ Source or Header  |  1994-03-19  |  5.6 KB  |  281 lines
  1. /************************************************************************
  2.  *                                    *
  3.  *                N O T I C E                *
  4.  *                                    *
  5.  *            Copyright Abandoned, 1987, Fred Fish        *
  6.  *                                    *
  7.  *    This previously copyrighted work has been placed into the    *
  8.  *    public domain by the author (Fred Fish) and may be freely used    *
  9.  *    for any purpose, private or commercial.  I would appreciate    *
  10.  *    it, as a courtesy, if this notice is left in all copies and    *
  11.  *    derivative works.  Thank you, and enjoy...            *
  12.  *                                    *
  13.  *    The author makes no warranty of any kind with respect to this    *
  14.  *    product and explicitly disclaims any implied warranties of    *
  15.  *    merchantability or fitness for any particular purpose.        *
  16.  *                                    *
  17.  ************************************************************************
  18.  */
  19.  
  20.  
  21. /*
  22.  *  FUNCTION
  23.  *
  24.  *    acos   double precision arc cosine
  25.  *
  26.  *  KEY WORDS
  27.  *
  28.  *    acos
  29.  *    machine independent routines
  30.  *    trigonometric functions
  31.  *    math libraries
  32.  *
  33.  *  DESCRIPTION
  34.  *
  35.  *    Returns double precision arc cosine of double precision
  36.  *    floating point argument.
  37.  *
  38.  *  USAGE
  39.  *
  40.  *    double acos (x)
  41.  *    double x;
  42.  *
  43.  *  REFERENCES
  44.  *
  45.  *    Fortran IV-plus user's guide, Digital Equipment Corp. pp B-1.
  46.  *
  47.  *  RESTRICTIONS
  48.  *
  49.  *    The maximum relative error for the approximating polynomial
  50.  *    in atan is 10**(-16.82).  However, this assumes exact arithmetic
  51.  *    in the polynomial evaluation.  Additional rounding and
  52.  *    truncation errors may occur as the argument is reduced
  53.  *    to the range over which the polynomial approximation
  54.  *    is valid, and as the polynomial is evaluated using
  55.  *    finite-precision arithmetic.
  56.  *    
  57.  *  PROGRAMMER
  58.  *
  59.  *    Fred Fish
  60.  *
  61.  *  INTERNALS
  62.  *
  63.  *    Computes arccosine(x) from:
  64.  *
  65.  *        (1)    If x < -1.0  or x > +1.0 then call
  66.  *            matherr and return 0.0 by default.
  67.  *
  68.  *        (2)    If x = 0.0 then acos(x) = PI/2.
  69.  *
  70.  *        (3)    If x = 1.0 then acos(x) = 0.0
  71.  *
  72.  *        (4)    If x = -1.0 then acos(x) = PI.
  73.  *
  74.  *        (5)    If 0.0 < x < 1.0 then
  75.  *            acos(x) = atan(Y)
  76.  *            Y = sqrt[1-(x**2)] / x 
  77.  *
  78.  *        (4)    If -1.0 < x < 0.0 then
  79.  *            acos(x) = atan(Y) + PI
  80.  *            Y = sqrt[1-(x**2)] / x 
  81.  *
  82.  */
  83.  
  84. #include <stdio.h>
  85. #include <math.h>
  86. #include "pml.h"
  87.  
  88.  
  89. #if !defined (__M68881__) && !defined (sfp004)        /* mjr++    */
  90.  
  91. static char funcname[] = "acos";
  92.  
  93. double acos (x)
  94. double x;
  95. {
  96.     struct exception xcpt;
  97.     double y;
  98.     
  99.     if ( x > 1.0 || x < -1.0) {
  100.     xcpt.type = DOMAIN;
  101.     xcpt.name = funcname;
  102.     xcpt.arg1 = x;
  103.     if (!matherr (&xcpt)) {
  104.         fprintf (stderr, "%s: DOMAIN error\n", funcname);
  105.         errno = EDOM;
  106.         xcpt.retval = HUGE_VAL;    /* for now, should ne NaN    */
  107.     }
  108.     } else if (x == 0.0) {
  109.     xcpt.retval = HALFPI;
  110.     } else if (x == 1.0) {
  111.     xcpt.retval = 0.0;
  112.     } else if (x == -1.0) {
  113.     xcpt.retval = PI;
  114.     } else {
  115.     y = atan ( sqrt (1.0 - (x * x)) / x );
  116.     if (x > 0.0) {
  117.         xcpt.retval = y;
  118.     } else {
  119.         xcpt.retval = y + PI;
  120.     }
  121.     }
  122.     return (xcpt.retval);
  123. }
  124.  
  125. #endif    /* !__M68881 && !sfp004    */
  126. #ifdef    sfp004
  127.  
  128. __asm("
  129. comm =     -6
  130. resp =    -16
  131. zahl =      0
  132. ");    /* end asm    */
  133.  
  134. #endif    sfp004
  135. #if defined (__M68881__) || defined (sfp004)
  136.  
  137.     __asm(".text; .even");
  138.  
  139. #ifdef ERROR_CHECK
  140.  
  141. __asm("
  142. _Overflow:
  143.     .ascii \"OVERFLOW\\0\"
  144. _Domain:
  145.     .ascii \"DOMAIN\\0\"
  146. _Error_String:
  147.     .ascii \"acos: %s error\\n\\0\"
  148. .even
  149. ");
  150.  
  151. #endif ERROR_CHECK
  152.  
  153. #ifdef ERROR_CHECK
  154.     __asm("
  155. | m.ritzert 14.12.1991
  156. | ritzert@dfg.dbp.de
  157. |
  158. |    /* NAN  = {7fffffff,ffffffff}        */
  159. |    /* +Inf = {7ff00000,00000000}        */
  160. |    /* -Inf = {fff00000,00000000}        */
  161. |    /* MAX_D= {7fee42d1,30773b76}        */
  162. |    /* MIN_D= {ffee42d1,30773b76}        */
  163.  
  164. .even
  165. double_max:
  166.     .long    0x7fee42d1
  167.     .long    0x30273b76
  168. double_min:
  169.     .long    0xffee42d1
  170.     .long    0x30273b76
  171. NaN:
  172.     .long    0x7fffffff
  173.     .long    0xffffffff
  174. p_Inf:
  175.     .long    0x7ff00000
  176.     .long    0x00000000
  177. m_Inf:
  178.     .long    0xfff00000
  179.     .long    0x00000000
  180. ");
  181. #endif ERROR_CHECK
  182.  
  183. __asm("
  184. .even
  185. .globl _acos
  186. _acos:
  187.     ");    /* end asm    */
  188.  
  189. #endif    /* __M68881__ || sfp004    */
  190. #ifdef    __M68881__
  191.  
  192.     __asm("
  193.     facosd    a7@(4), fp0    | acos
  194.     fmoved    fp0,a7@-    | push result
  195.     moveml    a7@+,d0-d1    | return_value
  196.     ");    /* end asm    */
  197.  
  198. #endif    __M68881__
  199. #ifdef    sfp004
  200.     __asm("
  201.     lea    0xfffa50,a0
  202.     movew    #0x541c,a0@(comm)    | specify function
  203.     cmpiw    #0x8900,a0@(resp)    | check
  204.     movel    a7@(4),a0@        | load arg_hi
  205.     movel    a7@(8),a0@        | load arg_low
  206.     movew    #0x7400,a0@(comm)    | result to d0
  207.     .long    0x0c688900, 0xfff067f8    | wait
  208.     movel    a0@,d0
  209.     movel    a0@,d1
  210.     ");    /* end asm    */
  211.  
  212. #endif    sfp004
  213. #if defined (__M68881__) || defined (sfp004)
  214.  
  215. # ifdef ERROR_CHECK
  216.     __asm("
  217.     lea    double_max,a0    |
  218.     swap    d0        | exponent into lower word
  219.     cmpw    a0@(16),d0    | == NaN ?
  220.     beq    error_nan    |
  221.     cmpw    a0@(24),d0    | == + Infinity ?
  222.     beq    error_plus    |
  223.     cmpw    a0@(32),d0    | == - Infinity ?
  224.     beq    error_minus    |
  225.     swap    d0        | result ok,
  226.     rts            | restore d0
  227. ");
  228. #ifndef    __MSHORT__
  229. __asm("
  230. error_minus:
  231.     swap    d0
  232.     moveml    d0-d1,a7@-
  233.     movel    #63,_errno    | errno = ERANGE
  234.     pea    _Overflow    | for printf
  235.     bra    error_exit    |
  236. error_plus:
  237.     swap    d0
  238.     moveml    d0-d1,a7@-
  239.     movel    #63,_errno    | NAN => errno = EDOM
  240.     pea    _Overflow    | for printf
  241.     bra    error_exit    |
  242. error_nan:
  243.     moveml    a0@(24),d0-d1    | result = +inf
  244.     moveml    d0-d1,a7@-
  245.     movel    #62,_errno    | NAN => errno = EDOM
  246.     pea    _Domain        | for printf
  247. ");
  248. #else    __MSHORT__
  249. __asm("
  250. error_minus:
  251.     swap    d0
  252.     moveml    d0-d1,a7@-
  253.     movew    #63,_errno    | errno = ERANGE
  254.     pea    _Overflow    | for printf
  255.     bra    error_exit    |
  256. error_plus:
  257.     swap    d0
  258.     moveml    d0-d1,a7@-
  259.     movew    #63,_errno    | NAN => errno = EDOM
  260.     pea    _Overflow    | for printf
  261.     bra    error_exit    |
  262. error_nan:
  263.     moveml    a0@(24),d0-d1    | result = +inf
  264.     moveml    d0-d1,a7@-
  265.     movew    #62,_errno    | NAN => errno = EDOM
  266.     pea    _Domain        | for printf
  267. ");
  268. #endif    __MSHORT__
  269. __asm("
  270. error_exit:
  271.     pea    _Error_String    |
  272.     pea    __iob+52    |
  273.     jbsr    _fprintf    |
  274.     addl    #12,a7        |
  275.     moveml    a7@+,d0-d1
  276. ");
  277. # endif ERROR_CHECK
  278. __asm("rts");
  279.  
  280. #endif /* __M68881__ || sfp004    */
  281.